function [ DrM ] = getDrM( r, rn, M )
%GETDRM Summary of this function goes here
%   Detailed explanation goes here

down = rn^4;

s = 1 / rn;
A = M;
ds = -r / rn^3;
dA2 = [2*r(1)*(r(2)^2 + r(3)^2), -r(2)*(r(1)^2 - r(2)^2 - r(3)^2), -r(3)*(r(1)^2 - r(2)^2 - r(3)^2), -2*r(1)^2*r(2), -r(1)*(r(2)^2 - r(1)^2 - r(3)^2), -2*r(1)*r(2)*r(3), -2*r(1)^2*r(3), -2*r(1)*r(2)*r(3), -r(1)*(r(3)^2 - r(2)^2 - r(1)^2);
       -r(2)*(r(1)^2 - r(2)^2 - r(3)^2), -2*r(1)*r(2)^2, -2*r(1)*r(2)*r(3), -r(1)*(r(2)^2 - r(1)^2 - r(3)^2), 2*r(2)*(r(1)^2 + r(3)^2), -r(3)*(r(2)^2 - r(3)^2 - r(1)^2), -2*r(1)*r(2)*r(3), -2*r(2)^2*r(3), -r(2)*(r(3)^2 - r(2)^2 - r(1)^2);
       -r(3)*(r(1)^2 - r(2)^2 - r(3)^2), -2*r(1)*r(2)*r(3), -2*r(1)*r(3)^2, -2*r(1)*r(2)*r(3), -r(3)*(r(2)^2 - r(3)^2 - r(1)^2), -2*r(2)*r(3)^2, -r(1)*(r(3)^2 - r(1)^2 - r(2)^2), -r(2)*(r(3)^2 - r(1)^2 - r(2)^2), 2*r(3)*(r(1)^2 + r(2)^2)] / down;
dA = reshape(dA2, 3, 3, 3);

DrM = tensorProductGradient(s, ds, A, dA);

end

